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Abstract 

Given a lattice basis of n vectors in Z n , we propose an algorithm using 12n 3 + 0{n 2 ) floating 
point operations for checking whether the basis is LLL-reduced. If the basis is reduced then the 
algorithm will hopefully answer "yes" . If the basis is not reduced, or if the precision used is not 
sufficient with respect to n, and to the numerical properties of the basis, the algorithm will answer 
"failed". Hence a positive answer is a rigorous certificate. For implementing the certificate itself, 
we propose a floating point algorithm for computing (certified) error bounds for the entries of the R 
factor of the QR matrix factorization. This algorithm takes into account all possible approximation 
and rounding errors. 

The cost 12n 3 + 0(n 2 ) of the certificate is only six times more than the cost of numerical 
algorithms for computing the QR factorization itself, and the certificate may be implemented using 
matrix library routines only. We report experiments that show that for a reduced basis of adequate 
dimension and quality the certificate succeeds, and establish the effectiveness of the certificate. 
This effectiveness is applied for certifying the output of fastest existing floating point heuristics of 
LLL reduction, without slowing down the whole process. 

1 Introduction 

Our motivation is to develop a certificate for lattice basis reducedness that may be used 
in cooperation with — possibly non certified — numerical reduction heuristics such as those 
described in [3T1 Ch. II-3] and [20] . The two main constraints are speed and effectiveness. 
Indeed, the certificate has to be fast enough for not slowing down the whole process, and 
the answer should be relevant ("yes") on a large class of inputs such as those successfully 
treated by the heuristic. Hence our general concern is somehow the compromize between 
speed and proven accuracy. The certificate will be introduced later below. It relies on error 
bounds for the R factor of the QR factorization of a matrix that we discuss first. 

Bounding errors for the factor R. Let A be an n x n invertible integer matrix. The QR 
factorization (see for instance [101 Ch. 19]) of A is a factorization A = QR in which the factor 
R g M nxn is an upper triangular matrix, and the factor Q G M nxn is orthogonal (Q T Q=T). 
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We take the unique factorization such that the diagonal entries of R are positive. Let F 
denote a set of floating point numbers such that the arithmetic operations in F satisfy the 
IEEE 754 arithmetic standard [1J. Assume that an approximate floating point and upper 
triangular factor R e F nxn is given. In Section [6] we propose an algorithm for computing a 
componentwise error bound for \R — R\ using operations in F only. For a matrix A = {ctij), 
\A\ denotes Our error bound for \R — R\ is given by a matrix H e F nxn with positive 

entries such that (see Qj on page [7]): 

\R-R\<H\R\. (1) 

Since floating point numbers are rational numbers, when R and E are known, (pQ) provides 
a rigourous mathematical bound for the error with respect to the unkown matrix R. 

For understanding the behaviour of the error bounding algorithm better, we recall in Sec- 
tion[3]some existing numerical pertubation analyses for the QR factorization. The necessary 
background material may be found in Higham's book [TU] . Then in Sections H] and 0, we 
give the mathematical foundations of our approach. We focus on the componentwise bounds 
of [3l] that allow us to derive an algorithm based on the principles of verification (self- 
validating) methods. On the latter methods we refer to the rich surveys of Rump [25, 26J, 
see also the short discussion in Sectional As numerical experiments of Section l6\3l will demon- 
strate, the error bounding algorithm is effective in practice. Its cost is only 5 times more 
than a numerical QR factorization, we mean lOn 3 + 0(n 2 ) operations in F. For efficiency, 
the error bounds are themselves calculated using floating point operations, nevertheless, they 
take into account all possible numerical and rounding errors. The reducedness certificate will 
require 2n 3 + 0(n 2 ) additional operations. Most of the 12n 3 operations actually correspond 
to the evaluation of matrix expressions. An efficient implementation may thus rely on fast 
matrix routines such as the BLAS [8]. 

At a given precision, the error bounding algorithm provides relevant bounds for input 
matrices with appropriate numerical properties. In particular, the dimension and related 
condition numbers should be considered in relation with the precision (see Section 16. 3j) . 
However, the power of the verification approach [23, is to be effective on many inputs for 
which the numerical approach itself is effective — here the numerical QR factorization. For 
example, we report experiments using 64 bits floating point numbers, and R computed by 
the modified Gram-Schmidt orthogonalization (see [101 Alg. 19.12]). On integer matrices of 
dimension n = 1500 with condition number around 10 5 , we certify that the relative error on 
the entries of R has order as small as 10 -6 or 10~ 5 , with only 10~ 10 or 10~ 9 on the diagonal. 
We refer here to the diagonal entries since they play a key role for instance in the LLL Lovasz 
test (see (EJ). For large condition numbers (with respect to double precision), say 10 12 , and 
n = 200, the algorithm may typically certify relative errors in 10 _1 , and 10~ 4 on the diagonal. 

The LLL-reducedness certificate. The effectiveness of the error bound on \R — R\ 
allows us to address the second topic of the paper. To an n x n integer matrix A we as- 
sociate the Euclidean lattice £ generated by the columns (aj) of A (for definitions and on 
algorithmic aspects of lattices we refer for instance to [B]). From (a,-), the LLL algorithm 
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computes a reduced basis [12], where the reduction is defined via the Gram-Schmidt or- 
thogonalization of a\, a^-, ■ ■ ■ , a n G Z n . The Gram-Schmidt orthogonalization determines the 
associated orthogonal basis * G Q n by induction, together with factors //jj, using 

a* = aj — Yl l j~=i^ij a *j^ an d = ( a i, a *j)/\\ a j\\h 1 — J ' < Vectors ai,a 2 ,...,a„ are said 
proper for 77 > 1/2 if their Gram-Schmidt orthogonalization satisfies 

|/%| < ^, 1 < j < i < n. (2) 

In general one considers 77 = 1/2. The basis a±, a2, ■ ■ ■ , a n of £ is called LLL- reduced with 
factors 5 and 77 if the vectors are proper, and if they satisfy the Lovasz conditions: 

(<*-/4i,*)IKIl2< \\ a * + i\\l 1 <<<»-!, (3) 

with 1/4 < <5 < 1 and 1/2 < 77 < \f5. If A = QR is the QR factorization of A then we have 

f IKH2 = r«, 1 < z < n, 

\ /iij = ^'i/^i, 1 < j < i < n. 

We see from (T4]) that if an approximation R of R with error bounds on its entries are known, 
then (depending on the quality of the bounds) it may be possible to check whether (T5]) and ([3]) 
are satisfied. All the above draws the reducedness certificate that we propose in Section [71 
We also fix a set F of floating point numbers, and perform operations in F only. For certifying 
the reducedness of the column basis associated to A the certificate works in three steps: 

I: Numerical computation of a R factor R such that A w QR; 

II: Certified computation of F G F nxn such that \R - R\ < F (see (P); 

III: Certified check of properness ([2]) and Lovasz conditions ([3]). 

Following the principles of verification algorithms [25], Step I is purely approximation, 
and we propose an implementation of Steps II and ill that is independent of the factorization 
algorithm used for computing R. For taking into account all possible numerical and rounding 
errors, Steps II and ill use certified computing techniques (see Section 16. ip . We rely on 
the fact that the arithmetic operations +, — , x, in F are according to the IEEE 754 

standard. We especially use explicit changes of rounding mode for certified bounds. 

Verification algorithms are a powerful alternative between numerical and computer alge- 
bra algorithms, they somehow illustrate the boundary between the two fields. The reduced- 
ness certificate we propose illustrates a cooperation of purely numerical computation with a 
certified approach based on the IEEE 754 standard, in order to provide a computer algebra 
answer. Our progress in linear algebra is in the line of previous works on error bounds for 
linear systems [2_2J EH EE] , on certifying the sign of the determinant [221 E] > on verifying 
positive definiteness [27], or on eigenvalues [TBI, 124] . Our contribution is to establish the 
effectiveness of componentwise bounds for a whole matrix, propose a corresponding certified 
algorithm using fast verification techniques, and derive and test with experiments a certifi- 
cate for the LLL reducedness application. 
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Absolute value and matrix norms. We already considered above the absolute value of a 
matrix A = (a y -) denned by (|ay|). We write \A\ < \B\ if \a%j\ < \bij\- It is possible to check 
that if A = BC then \A\ < \B\\C\. We will use several matrix norms (see [TUl Ch. 6]) such as 
the Frobenius norm || ■ \\p or the 2- norm || ■ H2. We will also especially use the infinity norm 
|| ■ ||oo = maxi<j<„ YJj=i \ a ij\- For A = BC we have ||A||oo < || J B|| o||C'|| 00 , and if h = WAW^ 
then \A\ < H with hy = h. 

Condition numbers. For a nonsingular matrix A, the matrix condition number is defined 
by K p (A) = ||A||p||A _1 || p with p = 2, F or 00 [TOl Th. 6.4]. With the infinity norm we will 
also use the Bauer-Skeel condition number cond(A) = |||j4 _1 ||>1||| C!C i < ftoo(^4) PHI §7.2]. 

2 Error bounds computation and verification algorithms 

In linear algebra, few things are known about the complexity of computing certified and 
effective error bounds. The problem is somewhere between the one of computing approximate 
solutions, and the one of computing multi-precision or exact solutions. A main result in [7] 
shows that the problem of computing a certified estimation of ||^4 _1 || (for a consistent matrix 
norm) is as difficult as testing whether the product of two matrices is zero. Hence if we 
consider 0(n 3 ) operations for multiplying two matrices of dimension n, a deterministic error 
bound — based on a condition number bound — would cost 0(n 3 ). The use of randomization 
may lead to error estimations in 0(n 2 ) operations, we refer to [TO], Chap. 15] and references 
therein, and to the fact that the matrix product could be verified in 0(n 2 ) operations [9]. 
We did not investigate the randomization possibilities yet. 

Verification methods have been developped in [231 El] for computing certified error bounds 
for linear system solution. In [21] the error bound (normwise) is computed in twice the time 
of numerical Gaussian elimination. In the same spirit, a verification approach using 0(n 3 ) 
floating point operations is proposed in [22] for the sign of the determinant (see [11] for 
a survey on this topic). Note that computing the sign of the determinant corresponds to 
knowing the determinant with a relative error less than 1. Our error bounding algorithm 
for R will also use 0(n 3 ) floating point operations. The verification approach [25], |26] gives 
an effective alternative to interval arithmetic whose exponential overestimation of the error 
would not be appropriate for our problem [26], §10.7]. The general strategy for calculating 
an error bound is first to establish a result whose assertion is a mathematical expression for 
the bound (see Theorem 14.21) . then design an algorithm that verifies the assumptions for the 
latter assertion, and computes a certified evaluation of the bound (see Section [5]). 

3 Perturbation analyses and bounds for the QR factorization 

A finite precision computation of the QR factorization of A leads to an approximate factor R. 
The errors in R with respect to R are called the forward errors (absolute or relative). The 
matrix R is not the factor of the QR factorization of A, however, it is seen as the QR factor 
of a perturbed matrix A = A + E, where E is called the backward error. The choice of 
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A is non unique, and one refers for instance for the smallest error norm. The link between 
backward and forward error is made using the condition number of the problem, hence for us 
the condition number for the problem of computing R. The (relative) condition number of 
the problem — under some class of perturbations — measures the relative change in the output 
for a relative change in the input. In this context, a useful tool for estimating the accuracy 
of the solution to a problem, is the rule of thumb [TOj p. 9] : 

forward error < condition number x backward error. (5) 

We survey below some more precise instantiations of (jSJ) for the QR factorization. Known 
results are, in general, approximate inequalities (first order results), but could be extended 
for giving strict bounds on the forward error. The rule of thumb therefore gives a first possible 
direction for deriving an error bounding algorithm for \R — R\ (the forward absolute error). 
However, most of corresponding bounds rely on matrix norms, and may thus overestimate 
the actual componentwise error in most cases. 

We will investigate an alternative direction in Section HI Rather than on the rule of 
thumb, our error bounding algorithm will be based on the componentwise bounds of Sun [34|. 
This will lead to an algorithm that seems to be naturally more effective than a matrix norm 
approach for our problem. Another advantage of using Sun's results is to remain in the spirit 
of the verification methods. In particular, we will see that the error bounding algorithm is 
oblivious of the algorithm that is used for computing the approximate factor R. Our bound 
computation may be appended to any numerical QR algorithm, and does not rely on back- 
ward error bounds that would be have been needed for using §5$). An approximate Q in not 
orthogonal in general, the backward error problem is to know for which matrix A close to 
A, there exists an orthogonal Q such that A = QR1 Backward error bounds are known 
for specific QR algorithms such as Householder or Gram-Schmidt ones (see Theorems 19.4 
and 19.13 in [10]), but may not be available in the general case. We will circumvent the need 
of the backward error in Section H] using the correspondence between the QR factorization 
of A, and the Cholesky factorization R T R of A T A. 

Sensitivity of the QR factorization. The condition number of the problem of computing 
R (the "rate of change" of R) in the QR factorization may be defined theoretically for given 
classes of perturbations, but it is non trivial to derive expressions of the condition number 
that can be used in practice. Nevertheless, various formulae are proposed in the literature 
providing quantities that can be thought as a condition number for R, we refer for instance 
to [I]. These quantities may be very effective in practice in a matrix norm setting. 

Let A = QR and A « A + E = QR be QR factorizations. As already noticed, for a 
floating point factorization A « QR, in general we have Q ^ Q since Q is not orthogonal. 
Let R = R + F. For a sufficiently small backward error E, consider the normwise relative 
error e = ||-E||i?/||y4|| 2 = \\A — A\\f/\\A\\2. Then Sun's [33j Rem. 3.5] perturbation bounds 
(see also [32]) give 

\\R - R\\ F /\\R\\ 2 < V2K 2 (A)e + 0(e 2 ). (6) 
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An improved bound is given by Zha [Theorem 2.1] |35j (see also [H §5] and [TQl §19.9]) under 
a componentwise model of perturbation that we simplify here. Let \A — A\ = \E\ = e\A\, 
then for sufficiently small e we have: 

\\R - i?||oo/||i?||oo < c n cond( J R- 1 )e + 0(e 2 ) (7) 

where c n is a constant depending on n. Hence the Bauer-Skeel condition number of R~ l can 
be considered as a condition number for the problem of calculating R. This indicates that one 
may potentially loose significant digits (in the result) linearly with respect to the increase of 
logcond(i? _1 ). This typical behaviour is illustrated by Figure 3.1 where we have computed 
QR factorizations of random matrices (of randsvd type [TUI Ch. 28]). The algorithm used is 
the Modified Gram-Schmidt algorithm [TDJ Algo. 19.12]. 
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Figure 3.1: Maximum relative diagonal error in R (Modified Gram-Schmidt algorithm) 
with respect to cond(ii ) for random matrices A (n = 200). 

Identities ([6]) and (J7J) provide first order estimations of the errors. They are essential for 
an idea of the normwise loss of accuracy. Nevertheless, the loss of accuracy on individual 
entries (needed for the reducedness certificate) may not be deduced from these identities. 
Consider for instance the case of Figure 3.1 where the ratios of the may be as large as 10 11 . 
The normwise bound of (JTj), that involves the max row sum ||i2||oo, cannot provide relevant 
informations for every |f»j — Tij\. Note also that the loss of accuracy would certainly be 
amplified by the implementation of the error estimation itself in finite precision (Figure 3.1 is 
a mathematical representation of the error). Normwise bounds much sharper than ([6]) and (JTJ) 
may be found, especially in [Sill], it remains to know how well the corresponding proposed 
estimations approximate the true condition number [H §10]. It would also be interesting to 
investigate how the new techniques of [SI H] could lead to practical componentwise bounds. 
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4 Strict componentwise bounds for the R factor 

We now present the mathematical view and justification of the error bounding algorithm 
of Section El Given A G ~R nxn invertible, and an upper triangular matrix R 6 M nx?1 , the 
problem is to bound \R — R\ where R is the unknown QR factor of A. In practice we will 
have A, Re ¥ nxn . 
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4.1 QR and Cholesky factorization 

The strict componentwise analysis of Sun [3H §4] for QR uses the matrix A such that 
A = QR is a QR factorization. Note that, because of the loss of orthogonality, if Q is a 
numerical approximation of Q then A is not in general the matrix QR. Informations on A 
may be available by taking into account the algorithm that has produced R. We refer for 
instance to [H Eq. (5.8)] and pTJl §19.9] where properties of Householder transformations are 
used for bounding the backward error. This is not sufficient for our problem since we are 
given only A and R, and since one of our goal is to be oblivious of the method used for R. 

For not relying on A, we propose to rather resort to Sun's study of the Cholesky factor- 
ization [3H §4]. If B G ~R nxn is symmetric positive definite, then there is a unique upper 
triangular R G ~R nxn with positive diagonal entries, such that B = R T R. This factorization 
is called the Cholesky factorization [TO, Th. 10.1]. It holds that A = QR is a QR factor- 
ization if and only if B = A T A = R T R is a Cholesky factorization. It may not be a good 
idea to use the Cholesky factorization for computing R numerically. The condition number 
of the problem may indeed increase too much, especially k,2{A t A) = (/t 2 (A)) 2 . For avoiding 
this drawback, our point is to implement the reduceness certificate of Section [7] using QR for 
computing R, and to use the Cholesky point of view only for computing the error bound. 

4.2 The bound on \R - R\ 

For a matrix A G IR nxn , the spectral radius p(A) is the maximum of the eigenvalue modules. 
We denote by triu(A) the upper triangular part of A, we mean that triu(A) = (tij) with 
tij = ciij if i < j, and tij = otherwise. The following Theorem is [MJ Th. 2.1]. 

Theorem 4.1. For B,B G IR nxn symmetric positive definite matrices, let R and R be the 
Cholesky factors of B and B. Let E = B — B , and 

G=\R~ T ER~ 1 \. (8) 

Then if p{G) < 1 we have 

\R-R\ < triu(G?(J — G?) _1 )|^|. (9) 

Inequality (Q is what we announced with ([TJ). Let us apply Theorem 14.11 with B = A T A 
and B = A T A. Using A = QR and Q T Q = I, we get from flS]): 

G = \Rr T ER-J\ = \R~J(B - B)R~ l \ = \Rr T (^A — A T A)R~ 1 \ ^ 

= \R~ T A T AR- 1 - R~ T A T AR~ 1 \ = \Q T Q - R~ T A T AR~ l \ = \R~ T A T AR~ l - I\. 

Going back to the R factor of the QR factorization we then have the following corollary to 
Theorem 14.11 

Theorem 4.2. For A G ~R nxn an invertible matrix, let R be the QR factor of A. Let 
R g ~R nxn be upper triangular and invertible, and 

G = \RT T A T ART 1 - I\. (10) 
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Then if 
we have 



P(G) < 1, 



\R-R\ < triu(G(J - G) )\R 



(11) 
(12) 



Proof. Since R is invertible, B = R T R is positive definite, the same holds for B = A T A. By 
construction R and R are the Cholesky factors of B and B. It suffices to apply Theorem 14.11 
for concluding. □ 

Few things are known about the (mathematical) quality of Bound ( fT2l) over IR. Further- 
more, both additional method and arithmetic errors will be introduced for the finite precision 
evaluation of the bound. Additional method errors will be introduced especially for calcu- 
lating certified bounds for R~ l and triu(G(J — G)~ x ) (see Section [5]). Additional arithmetic 
errors will be introduced by the finite precision itself. All together we produce an error 
bounding algorithm that is not fully analyzed, the experiments of Section 16.31 will however 
give a precise idea of its practical behaviour and effectiveness. For illustrating Bound f[T2l 
over IR, let us consider some examples that show that Theorem 14.21 leads to accurate bounds. 
The calculations have been done in Maple [16], either exactly or with high precision, then 
rounded for the presentation. Let H = triu(G(J — G)~ r ) such that (fT2l) is \R — R\ < H\R\. 

On the matrices used for Figure 3.1 (randsvd, n = 200), with R computed using 64 bits 
floating point numbers via the Modified Gram-Schmidt algorithm, we typically get the follow- 
ing. For A with cond(i? _1 ) ~ 10 5 , the infinity norm of the error matrix is ||-ff||oo ~ 2 x 10~ 9 . 
This leads to the knowledge that R approximates R with (relative) accuracy ~ 10~ 10 . The 
accuracy of R is about 10~ 13 for the diagonal entries, and the diagonal error estimation is only 
in a factor of 2 from the true diagonal error. If cond(_R _1 ) ~ 4 x 10 13 then \\H ^3x 10 -3 , 
and R is known with accuracy about 10 -2 (2 x 10~ 5 on the diagonal). The ratio between 
the estimation and the true error is less than 4 on the diagonal. Again, we will certainly 
loose accuracy with our finite precision implementation, but keep a very satisfying overall 
behaviour. Consider also the matrix quoted from [U Eq. 5.4]: 



Ai 



1 i-io- 10 

1 1 + 10~ 10 



with cond(-R *) ~ 2 x 10 10 . We compute the matrix R in Matlab |14j . and obtain over 
the error bound: 



\R-R\ 



9.7 x 10~ 17 -1.3 x 10~ 16 
3.7 x 10- 17 



< 



3.5 x 10~ 12 




3.5 x 10~ 12 
7.4 x 10~ 17 



(13) 



The matrix R is known with (relative) accuracy about 2.5 x 10~ 12 on the first row, and 
5.25 x 10~ 7 for r 22 . On the first row the error is overestimated by a factor about 3.6 x 
10 4 . Notwithstanding the fact that the accuracy of the bound produced by Theorem 14.11 is 



S 



penalized by the particular form of the matrix, the estimation of the accuracy of R remains 
very good. Now let A be the random 3x3 integer matrix 





-60 


28 


51 


A 2 = 


-24 


-35 


-89 




37 


51 


-23 



We look at Bound (TT2"j) when perturbing only the second row of the exact R and get: 



\R-R\ 




0.0071 -0.0052 




< 




0.014204 






0.023087 
2.9 x 10~ 6 



(14) 



The estimator computes the errors very well on the first and the second row. We think that 
the dummy error estimated for is a repercussion of the perturbation of row two. In next 
section we review the different quantities that are involved in Theorem 14.21 with the aim of 
looking at first implementation aspects. 



5 Toward an implementation 

Theorem 14. 21 is the foundation of our error bounding algorithm. It involves several quantities 
that need further study before deriving an implementation in Section [61 We decompose the 
computation of the bound on \R — R\ into four principal tasks. We need to: 1) check 
that R is invertible; 2) compute a bound on G; 3) check that p(G) < 1; and 4) bound 
H = triu(G(J — We recall that at this point, only A and R are known. 



5.1 Invertibility check of R 

For dealing with R~ l in a certified way, which is clearly a non trivial question in finite 
precision, we use the verification solution of Oishi and Rump [21]. We compute a purely 
numerical approximate inverse V ~ R~ l (by numerical triangular inversion). Then we know 
from [21] that, if 

\\RV - IWoo < 1, (15) 

then R is invertible. 



5.2 Bounding G 

For bounding G, and dealing with the unknown inverse of R, we are also inspired by [21J, 
and introduce W = RV (~ I). We have 

G = \R- T A T AR-^ - I\ 

= \(W- T W T )R- T A T AR- 1 (WW- 1 ) - (W- T W T )(WW- l )\ 

= \W- T (V T A T AV - W T W)W- V )\ < \W~ T \ ■ \V T A T AV - W T W\ ■ IPF" 1 ). 
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In the inequality above, if R is close to R and V is close to R 1 , then both V T A T AV and 
W T W are close to identity. Hence it is natural to pursue with: 

G < \W~ T \ ■ \V T A T AV - I + I -W T W\ ■ \W~ l \ 
< \W~ T \ ■ \(V T A T AV -I)- (W T W - I)\ ■ \W^\ 

which gives 

G < \W- T \ ■ (\(V T A T AV - I)\ + \(W T W - 7)|) • \W- l \. (16) 

We will use ffl6|) for computing a certified bound for G. The products involving A, R, 
V, and W = RV will be bounded directly by interval techniques. It remains to bound 
| 1 1 . We expect W to be close to I, and may use a specific approximation. We have 
\W~ l \ = |(I - (J - W))~ l \ (see |2H Intro.]). Then, when R is invertible, 

\W~ l \ =\I +(I -W) + (I -W) 2 + ...\ 

= \2I — W + (I — Wf{I + (I — W) + (I — w) 2 + ...)| 

< \2I -W\ + \(I - Wf\ ■ \I + (/ - W) + (/ - Wf + . . . | 

< \2I -W\ + M(\\I - W\\ 2 J{1 - \\I - WWoo)) 

where Ai(x) for x G M denotes the matrix whose all entries are equal to x. Here we have 
used the fact that the entries of \I - W\ 2 ■ \I + (I — W) + (I - W) 2 + . . . | are bounded by 
the infinity norm. Since W is triangular, it follows that 

IW- 1 ] < \2I -W\+ Jjlz^l- . triu(l n • l T n ) (17) 

where 1„ is the column vector with all entries equal to 1. Note that the invertibility check fTToT) 
ensures that 1 — \\I — WW^ > 0. The absolute value (H^ -1 ! could have been bounded directly 
using 1/(1 — || I — Wlloo), but introducing the infinity norm only in the second order terms 
leads to a much better bound in our experiments. 

The matrix manipulations we have done for obtaining ( fl6|) and ( TlTl) follow some keys to 
the design of verification methods. We especially refer to [261 P- 211] where the introduction of 
small factors is recommended. We have introduced the matrices V T A T AV — I and W T W — I 
whose absolute bounds are expected to be small when R rs R and W ~ I. On the other 
hand, in (IT7j) . \2I — W\ is expected to be close to 7, and remaining terms are second order 
terms (see also the analysis for a in [2H §5]). 

5.3 Bounding the spectral radius of G 

For any consistent matrix norm we have p(A) < \\A\\. With the above bound on G, we will 
simply test whether 

||G|| 00 <1 (18) 

for asserting that p(G) < 1 in Theorem l4.21 This test corresponds to the Gershgorin disks. It 
could certainly be sharpened in future versions of the certificate, see for instance the Cassini 
ovals in [3], or the iterative estimation in [27] . 
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5.4 Bounding \R — R\ 

Once a bound on G is known it remains to bound H = triu(G(J — G)^ 1 ). We have 

G{I - G)- 1 = G + G 2 + G 3 + ... = G + G 2 (I + G + G 2 + . . .) 
and 2 

triu(G(/ - G)- 1 ) < triu(G) + triu ' X " ' X ™) " ( 19 ) 

Since G is expected to be small, H = triu(G(J — G)^ 1 ) is expected to be close to triu(G). 
Note that using the spectral radius check ( |T8l) ensures that 1 — HGHoo > 0. 



6 Error bounding algorithm for the QR factor R 

Let F be a set of floating point numbers such that the arithmetic operations in F satisfy 
the IEEE 754 standard. A and R are now matrices in F nxn . Since (finite) floating point 
numbers are rational numbers, A and R can be seen as rational matrices. Let R G W nxn be 
the unknown QR factor of A (in general, the entries of R are not in F) . We carry the approach 
of Section [5] over to the floating point case for computing a floating point matrix H such that 
\R — R\ < H\R\. The error matrix H provided by Theorem 14.21 can be computed modulo 
the two checks (IT51) and (fTSl) . and using the inequalities (ITo]) . (117)1 . and (IT5|) . These checks 
and inequalities only involve matrix multiplications, additions, subtractions, and divisions 
by a scalar. After explaining the basic techniques we use for computing certified bounds 
in floating point arithmetic, we present the error bounding algorithm and demonstrate its 
effectiveness on various examples. 



6.1 Certified bounds for floating point matrix expressions 

We denote by fl(x) the value of an arithmetic expression x computed by floating point 
arithmetic in F. For instance, for a, b G F, fl(a + b x c) denotes the result in F with 
the addition and the mutiplication performed in floating point arithmetic. In the text, an 
arithmetic expression on floating point numbers denotes the exact value in K.. For instance 
a + b G R is the result of the addition in R. The abolute value, the max, and the negation 
are exact operations: for a, b G F, fl(|a|) = \a\, fl(max{a, b}) = max{a,6}, fl(— a) = —a. 

Thanks to the IEEE 754 standard, we can use the possibility of changing the rounding 
mode for computing certified bounds. We essentially follow Rump's approach for imple- 
menting verified matrix operations [261 . and Oishi and Rump [21]. We use the statements 
"setround(down)" and "setround(up)"Q All operations after a statement "setround(down)" 
or "setround(up)" are rounded downwards or upwards, respectively, until the next call to 
setround. For two floating point numbers a and b, a bound r on \a op b\ for op G {+, — , x , -=-} 
may be computed as follows. The program 

setround(down); r = H(aopb) , . 

setround(up); r = fl(a op b); r — max{|r|, |r|} 

* f esetround(FE_DDWNWARD) and f esetround(FEJJPWARD) in C language. 
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leads to r and r such that r < a op b <r, and to r G F such that \a op b\ < r, for any a and 
b, and any op. The IEEE standard ensures that r and r are the best possible bounds in F. 
This may be extended to the matrix operation A x B — C with A,B,C G F nxn . If A x B is 
implemented using only additions and multiplications, then the program 

setround(down); R = &(A x B — C) , . 

setround(up); S=fl(AxB-C); = max{|^|, ^ 

where the maximum is taken componentwise, provides R< Ax B — C < R, and R G W nxn 
such that | A x B — C\ < R. For bounding more general matrix expressions we will use a 
midpoint-radius matrix representation (we refer to [251 §10.9]). Assume that M and iV are 
two matrices known to be in intervals [M, M] and [N_, N] , respectively. The intervals are for 
instance obtained by a computation of type (j2"T|) . Then the program [261 Fig- 10.22]: 

setround(up); m M = fl((M - M)/2); r A/ = fl(m A/ - M) 

mjv = fl((JV-i\0/2); r^^fl^-iV) 
setround(down); R = &(m.M x nijy — /) (22) 
setround(up); i? = fl(m M x m^r — I) 

R = ft (max{\R\, \R\} + |m M | x r^ + r M x (|mjv| + r N )) 

computes R such that \M x N — I\ < R. Both (j2*T|) and ff22l allow to use fast matrix routines 
such as the BLAS ones (see the general discussion in [26l §10.9]) The number of operations 
in F needed is 2 and 4 matrix products, respectively. 

Other matrix operations that we will perform are additions, products, and divisions by 
scalars for matrices with positive entries (absolute values essentially). We also compute in- 
finity norms. With no subtraction involved, certified bounds can be computed using directed 
rounding. From (120]) . upper bounds for these computations are obtained by evaluating the 
floating point expressions after a "setround(up)" statement. For upper bounds on divisions 
by a floating point number 1 — g, we first compute upper bounds for —(g — 1) and l/(g — 1). 

Other approaches for certified matrix computations could be considered. We refer to 
Rump [26] for a general discussion on this topic, and for the efficiency of the approach 
chosen here. 

6.2 Computing an error bound 

For A and R in F nxn , R upper triangular, we follow Section for computing a floating point 
matrix H such that \R — R\ < H\R\. All operations are done in the given floating point 
number set F. For simplifying the presentation we often forget the costs in 0(n 2 ). 

The first step is the computation of V ~ R" 1 . Such a triangular matrix inversion is done 
in n 3 /3 operations [101 Ch. 14]. We then compute W_ and W for W = RV by two triangular 
matrix products, this is done in 2n 3 /3 operations. This dominates the cost for checking that 
R is invertible by bounding \ W — 1\ using (I2T1) . and by the infinity norm test (|15|) . Re-using 
W and W, a bound on l^ -1 ] is then computed using (1T71) in 0{n 2 ) operations. The latter 
uses (121]) for \ W — 2I\, and computes a bound with positive matrices using directed upwards 
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rounding. We now look at bounding G using (|T6|) . Since G is symmetric we restrict ourselves 
to counting the operations for calculating the upper triangular part. With W G [W, W) one 
can bound \W T W — I\ using (1221) in four matrix products. Since W is upper triangular, 
and W T is lower triangular, the bound is obtained in 4n 3 /3 operations. We then use ([21]) 
and (1221) for computing an interval for AV in 2n 3 operations (two dense x triangular matrix 
products), and for bounding \ V T A T AV T — 1\ in 4n 3 operations (four dense products resulting 
in a symmetric matrix). A bound on G is deduced by operations on matrices with positive 
entries in 4n 3 /3 operations. The latter is essentially two dense x triangular matrix products 
with a symmetric result. Once a bound on G is known, testing its spectral radius by (fl8l) 
costs 0(n 2 ) operations. G has positive entries, a bound on the error matrix H can then be 
computed by directed towards rounding using (fT9l) also in 0{n 2 ) operations. 

We summarize this analysis, and take into account the final matrix product H\R\ in the 
following result. 

Theorem 6.1. Let A G ^" nxn j and R G W nxn upper triangular be given. The error bounding 
algorithm computes a matrix F G W nxn such that \R — R\ < F, where R is the unknown QR 
factor of A, in 10n 3 + 0(n 2 ) floating point operations. 

A QR factorization typically costs 2n 3 + 0(n 2 ) (Gram-Schmidt or Householder ap- 
proaches) or 3n 3 + 0(n 2 ) (using Givens rotations). Hence we are able to compute a cer- 
tified error bound \R — R\ at the cost of only five approximate factorizations. We have 
implemented the algorithm in C language. The error bounding program takes in input two 
floating point matrices A and R and always returns a matrix F. The entries of F are finite 
(positive) floating numbers if the program is able to certify that R is invertible, that the 
spectral radius of G is less than one, and if no overflow is produced. Otherwise, the entries 
of F may be equal to infinity. 



6.3 Computational results 



The results we present here correspond to the application of Theorer d6 . 1 1 wit h 64 bits floating 
point numbers. In this section and in Section [7| the condition numbers and the "true errors" 
have been computed with high precision using Mpfr [18J. For several types of matrices, we 
study the behaviour of the certified error bound by looking at its value and its accuracy 
(with respect to the true error), especially when the dimension and the condition number 
increase. We mainly focus on the exponent k such that relative error is in 10~ fc , k expresses 
the number of significant decimal digits we certify for the entries of R. Let us first come 
back on the examples of Section HI On the matrix Ai, and R from Matlab, we compute the 
bound 

6.7 x KT 11 6.7 x 10" 11 



R-R \< 







5 x 10 



-16 



Comparing to (fT3"j) . we see that the finite precision estimator we propose is only slightly 
overestimating the best bound that could be obtained by the method. On the matrix A2, 
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and the corresponding perturbation of the exact R we get: 

8.8 x 1(T 6 9.52 x KT 6 1.96 x 1(T 6 
/? - R\ < 0.014207 0.023098 

1.16 x 10~ 5 



The "large" perturbation of the second row is detected very accurately. For next results, R 
is computed with the Modified Gram-Schmidt algorithm using 64 bits numbers as for the 
estimator. Our tests use ten matrix samples. 
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Figure 6.1: Certified ||-ff||oo for random matrices A with K2(A) ss 10 3 . 

We first illustrate the value of the certified bound with respect to the dimension. Figures 6.1 
and 6.2 are for random input matrices A (of rands vd type [10J Ch. 28]). 
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Figure 6.2: Certified maximum relative error on R for random matrices A 
such that K2(A) « 10 3 (y axe with logarithmic scale). 

We keep the condition number almost constant when the dimension increase. We draw the 
infinity norm of H such that \R — R\ < H\R\, and the certified maximum relative error on 
the diagonal of R, we mean maxj |r# — Tn\/\ru\. We see that ||i?||oo increases linearly with n. 
The loss of accuracy on the diagonal is approximately quadratic in n (we use a logarithmic 
scale for the y axe on Figure 6.2). Such small increase rates — that are typical of numerical 
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algorithm forward errors themselves — demonstrate a first aspect of the effectiveness of our 
finite precision bounds. The certified general maximum error max^ jr^- — rjj|/|r^-| increases 
faster. It typically grows from 10~ 7 to for the dimensions considered here. We need 
further investigation for a better understanding of the latter behaviour, especially of the 
influence of the product H\R\, and of the magnitudes in R. Note also that for the two latter 
figures, cond(.R -1 ) is sligthly growing, and the growth of the estimation depends on the true 
error itself. 

We discuss next the accuracy of the certified bound with respect to the exact error (not the 
quality of the QR algorithm itself). In addition to above randsvd matrices we also consider 
random integer matrices with entries of absolute values less than 1000. On these two types 
of matrices we obtain similar results. The condition numbers ^(A) are varying from about 
10 4 to 10 6 . On random integer matrices of dimension 1500, the maximum exact relative 
error on R has order 10~ 10 to 10~ 9 . We are able to certify this error by returning an error 
bound of order 10 -6 to 10 -5 . With respect to the dimension, we observe that the fast certified 
bound overestimates the componentwise error by a factor of order about 10 3 for n = 200 
to about 10 5 for n = 1500. Restricted to the diagonal entries, the overestimation goes from 
about 10 2 to less than 10 4 . This shows that even with condition numbers and dimensions 
that can be here quite large, we are able to certify at least four or five significant decimal 
digits for every entries of R, and at least 9 digits on the diagonal (where the error itself is 
much smaller in general). On matrices with small condition number (generated using Matlab 
gallery ( ' orthog' ) [TUJ Chapter 28]) the quality of the certified bound may be remarkably 
small and stable with respect to the dimension. For dimensions between 60 and 500, and 
cond(i?~ 1 ) ~ 3 (ftoo < 200), we most of the time obtain an overestimation between 15 and 
22 (and more than 12 certified significant decimal digits in R). 

We may now ask the question of the sensitivity of the quality of the certified error bound 
with respect to the condition number of the input matrix. We first report that the quality 
maybe be very good even for matrices with high condition number. For Figure 6.4 we 
use A = QAx £ F nxn . The matrices Q are random orthogonal from the Matlab gallery 
function [TUJ Chapter 28] . The matrices are Kahan upper triangular matrices with an = 
(sin^)* -1 , aij = — (sin^) 4 " 1 cos 6 for j > i, and 8 = 1.2. 



Dimension 


10 


20 


30 


40 


50 


60 


70 


Koci(A) 


10 2 


1.3 x 10 4 


1.1 x 10 6 


7.8 x 10 7 


4.8 x 10 9 


2.8 x 10 11 


1.5 x 10 13 


Bound / error 


45 
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281 
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152 


Certified digits in R 


14 


12 


10 


9 


7 


5 
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Figure 6.4: Ratio of the certified relative error bound and the true error (max) 
on Kahan matrices, and number of significant decimal digits certified in R. 

However, in general, the quality of the bound may depend on the condition number. Consider 
for instance the ratio of the certified relative error bound and the true error (max) for small 
matrices (n = 10). For a Chebyshev Vandermonde-like (nearly orthogonal, ~ 13), the 
ratio is about 11. We have a ratio about 14 for Toeplitz and symmetric positive definite 
matrices (/t^ ~ 700). On the Pascal matrix (k^ 8 x 10 9 ) we get a ratio about 25, and 
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about 1600 for the Hilbert matrix ( 10 13 ). Figure 6.5 is more general. The 

overestimation of certified error bound seem to increase quite slowly with the condition 
number. 

1600 j j j 1 

1 

§ eoo / 

o 

/ 

I 

1 400- ./ 



10 A 3 10*7 T0*11 Condition number 

Figure 6.5: Ratio of the certified relative error bound and the true error (max) 
with respect to k 00 (^4) on randsvd matrices of dimension n — 200. 

We see that the limits of our algorithm, we mean the conditions in which it is returning 
finite bounds, are clearly linked with the numerical properties of A. Let us give two examples 
for the impossibility to certify the spectral radius using (fl8|) . We return finite bounds for 
the error on every entries of R for the Pascal matrix of dimension 14 (k^ « 3.8 x 10 14 , 
||G||oo ~ 0.06). For n = 15 the algorithm produces infinity bounds. On random randsvd 
matrices of dimension 40, the algorithm is effective until ~ 3 x 10 14 with \\G\lao « 0.9. 
Note that in double precision, with relative rounding unit 2~ 53 (the backward error is larger 
in general), and for a relative forward error less than 1, the rule of thumb (jSJ) advocates for 
a condition number less than 10 16 . 

The certified bound is computed with finite precision, hence inherently, it overestimates 
the true error. However, for realistic dimensions and condition numbers (with respect to the 
precision), the overestimation is mastered. It follows that in general, many significant digits 
are certified in the approximate QR factor R. The latter is a key to the application of the 
fast bound to the reducedness certificate. 

7 A certificate for LLL reducedness 

To an n x n integer matrix A we associate the Euclidean lattice £ generated by the columns 
(a,j) of A. About lattices the reader may refer for instance to Cohen's book [5]. Since the 
seminal Lenstra-Lenstra-Lovasz algorithm [12J — whose range of application is exceptional — 
the lattice basis reduction problem receives much attention. In particular, floating-point 
variants that lead to very fast reduction approaches have been invented. See the work 
of Nguyen and Stehle [191 EI]) °f Schnorr [29], and references therein. Most of floating 
point variants lead to powerful heuristics, especially a la Schnorr- Euchner [30], that are 
implemented (often with improvements) in most of computer algebra and number theory 
systems. Our aim here is not to study the basis reduction itself. We focus on the reducedness. 
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Indeed, a fast heuristic may not certify that the output basis is reduced (still working very 
well), and it is worthwhile to study the problem of checking a posteriori whether a given 
basis is reduced or not. The notion of reduction we consider is the LLL reduction [12]. 

We propose here an algorithm that takes in input an invertible matrix A e Z nxn , and 
tests the LLL reducedness of the basis formed by the columns of A. In the Introduction 
we have seen that this consists in testing the two conditions ([2]) and Q. Let R be the QR 
factor of A. If the aj are proper, we mean 

\rij\/r iti <rj, 1 < i < j < n, (23) 

and if the Lovasz conditions 

\js- (r iti+1 /r it i) 2 r M < r i+lii+ i, 1 < i < n - 1, (24) 

are satisfied, then the basis a±, a 2 , ■ ■ ■ , a n of C is called LLL reduced with parameters (5, 77). 
The latter satisfy 1/4 < 5 < 1 and 1/2 < r] < V5. 

The principle of the algorithm is to compute an approximate R together with error 
bounds (using the floating point algorithm of Section E]) , then to test (1231 and (124|) . 

The entries of A are integers of arbitrary size (our implementation relies on Gmp |17j). 
Therefore the entries of A may not be represented exactly by elements in F. Nevertheless, 
for the computation of an approximate R we may take A by direct conversion to F. Since the 
error is very small and R will be an approximation anyway, this does not really influence the 
quality of subsequent computations. Then R is computed by the Modified Gram-Schmidt 
algorithm. Once R is known we apply Theorem 16.11 for computing a certified error bound. 
The only expression that has to be bounded with A involved is in (fl6|) . where the computation 
of AV using the program (12T|) is needed. The problem of conversion to F is solved here by 
rounding upwards and downwards during the conversion integer to floating point. We mean 
that we introduce a small interval such that A e L4_,A + ] with A_,A + e ]p™ x ™ ( see the 
certified techniques in Section I6TT1) . and we evaluate A_V and A + V in (I2"T|) . Therefore the 
error bound F G F nx " we compute by Theorem 16.11 is actually such that \R — R\ < F for R 
the QR factor of any A e L4_, A + }. 

Once F is known, for fixed i and j, we test (|23|) by resorting to the bounding techniques 
of Section 16.11 

setround(down); r/ = fl^); t L = fl((r M - f iA ) x 77) 
setround(up); t j = &(\r id \ + f id ) (25) 
test tj <t{! 

with temporary variables ti and tj. Recall that the diagonal entries of R are positive. 
Similarly, for a fixed i, we test f )24l) using: 

setround(up); U = fl(r^ + /*,*); 5 = fl(5); 
setround(down); t i+1 = fl(r i+M+1 - 

t = fl(((|r M+1 | - /i, i+ i)A<) 2 ) - S; t = -t; (26) 
setround(up); t = fl(v^ x U) 

test t < 
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with temporary variables t and ij. In practice, for minimizing the cost induced by the 
changes of rounding mode, loops are put between the setround instructions. In addition 
to the 10n 3 + 0(n 2 ) operations for computing F using Theorem 16.11 the reducedness test 
essentially requires 2n 3 + 0(n 2 ) operations for computing an approximate factor R. This 
gives the following. 

Theorem 7.1. Let A 6 Z nxn invertible and parameters (5, if) be given. The reducedness 
certificate certifies in 12n 3 + 0(n 2 ) floating point operations that the column lattice of A is 
LLL reduced with parameters (5,rf), or returns "failed". 

The reducedness is certified when the error bound computed for \R — R\ is finite, when no 
overflow or underflow occur during the test, and when the basis is reduced. The cost of the 
certificate is roughly the one of six floating point QR factorizations. Therefore in general, 
the reducedness test should be much faster than the reduction process itself, and may be 
appended to any reduction heuristic program. 

Let us now report some experiments. As previously in the paper all codes are run using 
64 bits floating point numbers. The effectiveness of the certificate essentially relies on the ef- 
fectiveness of the error bounding algorithm. We have manipulated lattices using Magma [13j, 
the LLL reduction implementation is based on the work of Nguyen and Stehle [191 El] • The 
first family of reduced bases — matrices A — we consider are obtained by the reduction of 
n x n random integer matrices. The bases are reduced for the classical LLL parameters 
[S,r]) = (3/4,1/2) in Figure 7.1, and (5, rf) = (0.99,0.5001) for a stronger reduction in 
Figure 7.2. 



Dimension 


40 


200 


500 


1000 


K oo {-A) 

tk—t = min 4 {^+i - t] in ([26]) 
Certified absolute error on \\a* k 2 


4.7 x 10 2 
18 

7.5 x 10~ 12 


2.4 x 10 4 
10 

3 x 10" 10 


1.8 x 10 5 
13 

1.5 x 10" 9 


9 x 10 5 
23 

1.2 x 10" 8 


Certified max,, fj,ij 
Max. certified relative error on |rjj| 


0.4997 
2.8 x 10~ n 


0.499994 
8.6 x lO- 9 


0.49991 
1.5 x 10" 7 


0.49999 
3 x 10~ 5 



Figure 7.1: Reducedness certificate output on (3/4, l/2)-reduced bases from random 
integer matrices with entries on 10 3 bits, max lay | < 1000. 



Since the numerical quality of the tested bases is good (/too {A) < 10 6 ), the reducedness 
certificate is highly efficient. We mean that the certified error is very small, and hence the 
tests are passed except in exceptional cases. Figures 7.1 and 7.2 for instance look at the 
smallest difference t). — t whose positiveness has to be certified in f)26p . The certificate has 
lots of room since the absolute errors on t and tk = \\alW2 are much smaller. Exceptional 
cases will rather occur when testing properness. Indeed, testing reducedness may be an ill- 
posed problem because of the possible equalities in fT23|) and fT24]) . An ill-posed case with say 
1] = 1/2, is for example a reduced basis with = 1/2 for some Therefore the algorithm 
will rather be used for certifying that a (5, ?y)-reduced basis is a (S — ex, r) + e 2 )-reduced basis 
for small ex, 62- The latter does really affect the relevant certified informations provided by 
the reduction. 
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Dimension 


40 


200 


500 


1000 


tk —t = minj{ti + i — t} in (|26[) 
Certified absolute error on ||a|||2 


4.8 x 10" 2 
9.4 x 10" 14 


7.7 x 10- 2 
6 x 10" 12 


5.3 x 10- 2 
4 x 10" 11 


7.3 x 10" 2 
2 x 10" 10 



Figure 7.2: Reducedness certificate output on (0.99, 0.501)-reduced bases from random 
integer matrices with entries on 10 bits, max \a,ij \ < 10. 



A second type of reduced bases on which we have run the certificate comes from the 
problem of computing a good floating point coefficient polynomial approximation to a func- 
tion [2]. We have considered reduced bases with parameters (3/4,1/2). These bases may 
have integer entries as large as 10 80 . The certificate has always succeeded. On a 18 x 18 
example, with k^A) ~4x 10 12 , the smallest difference t — % has been around 2.4 x 10 76 
with certified absolute error 1.95 x 10 62 . The maximum of the //^ has been certified to be 
less than 0.493. On an example with n — 31 and Koo(A) ~ 8 x 10 13 , we have certified an 
absolute error 3.2 x 10 53 for t — tk ~ 1.7 x 10 67 . On the latter example we have also checked 
that max < 0.4991, thanks to a maximum relative error \R — R\ certified to be less than 
0.2 (only 6 x 10~ 15 on the diagonal). 

The first main source of failure of the certificate is the failure of the error bounding 
algorithm when the precision is too small compared to the numerical quality of the tested 
basis. We have run the certificate on a third class of reduced bases. These bases are obtained 
by the reduction of "random" (knapsack type) lattice bases in the sense of [201 §3.4]. In the 
experiments reported here, the non reduced bases have random integers of 10 3 bits in the 
knapsack weight row. The reduced bases in input of the certificate (matrices A) are dense 
with integers as large as 10 for n = 75, and 10 20 for n = 300. We use the parameters 
(5,rj) = (3/4,1/2) and (6,rj) = (0.99,0.5001). The choice (5,r]) = (0.99,0.5001) produces 
better reduced bases as shown by in Figure 7.3 (for a same non reduced basis). Until 
dimension 175 the certificate is very likely to succeed since the maximum certified relative 
error is small. On several tenths of trials, the certificate never failed, with a certified max 
as close to 1/2 (with r) = 1/2) as 0.4999916. 



Dimension 


75 


100 


125 


150 


175 


(6,t,) = (3/4,1/2), Koo (A) 
tk - t = min t {t 1+ i - t} in ((26j) 
Max. certified relative error on (r^ | 


6 x 10 5 
1.3 x 10 37 
1.3 x 10~ 9 


5.2 x 10 B 
8.5 x 10 26 
3.4 x 10" 8 


2.3 x 10 s 
4.2 x 10 20 
2.2 x 10" 6 


1.3 x 10 10 
3 x 10 15 
2.1 x 10" 5 


2 x 10 11 

1.2 x 10 12 

6.3 x 10" 3 


(S, r,) = (0.99,0.5001), Koo (A) 
Max. certified relative error on |r,j| 


2.4 x 10 4 
5.1 x 10- 10 


4.6 x 10 5 
2.5 x 10" 9 


4 x 10 5 
3.9 x 10" 8 


4 x 10 7 
6 x 10~ 7 


9 x 10 8 
9.5 x 10" 6 



Figure 7.3: Reducedness certificate output on "random" reduced bases from knapsack problems, 
max |ajj| goes from 10 45 (n = 50) down to 10 25 (n = 175). 



Beyond dimension 175 with this type of reduced basis, the certificate starts to fail more 
often. On dimension 200 with a conditioning about 10 12 with (3/4,1/2), the error bound 
on the relative error approaches 1. The properness with 77 = 1/2 may become impossible 
to check, and ask for a certificate with rj = 1/2 + e, say rj = 0.5001. Note that the Lovasz 
test (1261) seems to fail later thanks to much better error bounds on the diagonal in general. On 
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dimension 300 for (3/4, 1/2) the quality of the reduced bases is too deteriorated (fi^ ~ 10 19 ), 
and the error bounding algorithm fails with the impossibility of having a small spectral radius 
in (15.31) . Nevertheless, on a typical example of dimension 300 with a (0.99,0.5001) reduced 
basis, the error bounding algorithm remains effective ( 10 13 , Halloo » 0.6). The 

certificate may not be able to certify the actual reducedness of the basis, for example with 
minjjij — t} ~ —4.12 x 10 8 , and a too big absolute error bound 4.42 x 10 8 . By changing the 
certificate parameters to (5 — 61,77 + e 2 ) = (0.985,0.515), the certificate succeeds again, and 
therefore is still able to certify a relevant information on the basis. 

The numerical limitations of the certificate are close to those identified in [20l Heuristic 4] 
for the reduction process itself. Indeed, on the knapsack bases, it is claimed in [20] that a 
precision n/4 + o(n) should suffice when using the floating point reduction of [19]. This 
means that n pd 200 is a barrier with a 53 bits precision (64 bits numbers). The eventuality 
of a link between both limitations deserves to be further investigated. 

8 Conclusions 

Between numerical approximation and computer algebra, we propose a certificate for an 
(exact) algebraic/geometric property — the LLL reducedness of a lattice basis. This work, 
based on the fast computation of certified error bounds, inherits from the verification methods 
approach. In particular, thanks to the IEEE arithmetic standard, the floating point errors 
do not put a curb on the objective of certification. They may rather be mastered and 
used for accelerating the programs. In error bound computation and property certification, 
the foreground of our study is to understand the compromize between the cost and the 
quality/effectiveness of bounds and certificates. In our case for instance, may we hope for 
an 0(n 2 ) effective certificate? Various computer arithmetics come in the background, where 
floating point computation, multi-precision, verification identities, midpoint-radius intervals, 
and exact computation are collaborative tools. 

We think that our study raises several directions that deserve further investigations. The 
error bounding problem for the R factor, and its finite precision implementation should be 
better understood and improved, ingredients such as diagonal scaling and other approximate 
QR factorizations may be introduced. The usefulness of taking into account the algorithm 
used for computing R should be studied (in a more restrictive verification approach). A 
more general question is to know whether reducedness could be certified without resorting 
to the QR factorization? 

To our knowledge, the minimum precision required for a proven LLL variants is 1.6n+o(n) 
with the L 2 algorithm of [13 [2D] (for 5 close to 1 and 77 close to 1/2). Our experiments 
show we may certify reducedness for dimensions much higher than this worst-case limit 
(n mm < 53/1.6 ~ 33). The certificate is therefore very effective for a use complementary 
to reduction heuristics in dimension greater than n max with double precision. Noticing the 
fact that the certificate is sensitive to the numerical properties of the input basis, it is worth 
studying its extensions to reduction algorithms and reducedness certificates with adaptative 
precision, and sensitive to the numerical quality. 
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